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1. Introduction 


We have completed an investigation of urbanization in the State College area and have 
begun work on two additional target sites: Chester Co., Pennsylvania, and Costa Rica. 
State College served as a test area for developing our methodology: the target area is 
small and lightly populated, but it is undergoing rapid urbanization. Chester Co., however, 
is much larger, and is riddled with urban centers, some of which are undergoing rapid 
development, -mostly at the expense of farmland; {ref.. Table 2, Section 2). Costa Rica is 
undergoing both a very high rate of urbanization around its principal city, San Jose’, and 
had been experiencing a very large overall rate of deforestation (about 4% a year) until a 
decline began a few years ago. 

The purpose of the investigation is to use a new method for deriving land surface 
parameters from a combination of thermal infrared and vegetation index measurements 
obtained by satellite (Landsat-TM and NOAA-AVHRR) and to integrate these parameters 
with more conventional data bases. Our intent is to show that it is possible not only to 
monitor land use changes over time by satellite with conventional descriptors (percent 
deve ied, agriculture, forest, water) but to express these same land use change in terms 
of phxsically based surface climate variables derived from the satellite measurements, — 
moisture availability (Mo) and fractional vegetation cover (Fr). By physically based we 
mean that the surface energy fluxes can be calculated from these two parameters, which 
are the dominant ones governing the surface energy budget in land surface models. 
Moreover, given a field of Mo and Fr and a specific urbanization or deforestation 
scenario, future changes in local microclimate can be made with the aid of a atmospheric 
mesoscale model. 

Underlying all this work is our basic hypothesis, which is that changes in land use, 
including deforestation, exert a profound influence on local microclimate whose effects 
may greatly exceed in importance those occurring on larger scales. We hope to establish 
for the first time a systematic and quantitative approach to classification of images using 



thermal infrared (and vegetation index) measurements by satellite, from which the urban 
planner, the meteorologist and the geographer can assess both past and future impact on 
the environment. 

2. Present Investigation 

a) A Review of Results from the State College study 

A major research effort in our research has been devoted to the study of urbanization 
using the so-called triangle method (Gillies and Carlson, 1995; Gillies et al., 1997). In the 
first phase, Owen (1995; see also Owen et al., 1996 and Carlson and Owen, 1996 
(attached)) completed an analysis of the State College, PA, area in which values of Fr and 
surface moisture availability. Mo,, were determined from AVHRR data over a 30x30 km 
domain centered on State College, PA. The object was to see if the migration of pixels 
within the Fr/Mo space (the triangle) were coherent and could be related to the nature and 
extent of urbanization. Eight AVHRR scenes were obtained (two per summer every other 
year) between 1986 and 1994. Additional Landsat TM scenes were obtained, from which 
percent urban development was calculated and classified as to forest, agriculture, 
developed and water. 

Results, summarized in Owen et al. (1996) and Carlson and Owen (1996), are that (1) 
AVHRR data can detect neighborhood-scale changes in urbanization, (2) moisture 
availability (Mo) seems to constitute an intrinsic descriptor of the surface fabric, in 
addition to being a fundamental land surface parameter which governs local microclimate, 
(3) changes in Fr and a scaled infrared surface radiant temperature (T* r ) (as well as 
evapotranspiration) accompany significant changes in urbanization; Mo, however, does 
not appear to change significantly with increasing urbanization, at least when agricultural 
land is developed; Fr decreased substantially from the summer of 1985 to the notoriously 
dry summer of 1988 (but recovered in 1991), whereas Mo hardly changed during that 
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period. (4) tree cover appears to be a dominant mitigating factor in the aforementioned 
changes in land surface parameters. Tree cover seems to affect Mo. 

b) Chester Co. 

Following our methods developed for State College, analyses have begun of images for 
Chester County, PA. Chester Co., located just west of Philadelphia, is about four times the 
area of the State College target area. Parts of the county are undergoing rapid 
urbanization. 

i 

Table 1 : List of Scenes Obtained for Chester Co. 

(AVHRR and Landsat Thematic Mapper) 


J AVHRR 

Landsat - TM 

| 7/24/85 

6/10/87 

| 7/5/86 

7/12/87 

7/16/88 

6/12/88 

7/12/89 

5/20/91 

8/2/90 

7/28/93 

9/1/90 

6/2/96 

7/10/91 


8/17/91 

8/22/92 

7/17/93 

8/24/94 

7/15/95 

1 7/16/95 


AVHRR scenes from 1986 to 1995 (Table 1) have been converted to temperature and 
radiance values, corrected for atmospheric attenuation and georeferenced; this has also 















been done for a fewer number of Landsat TM scenes covering roughly the same period. 
Initially, the procedures parallel those applied to the State College images. As in the case 
of State College, summer scenes are used because they are subject to the least variability 
due to changing meteorological conditions, and because summer scenes permit a better 
distinction between forest, agriculture and developed land than spring or winter scenes. 
From these image data, we determine Mo and Fr and then relate these parameters to the 
conventional land use categories as they evolve over the previous decade. 

All TM images have been classified according to several basic land surface types. 

Although the latter operation involves designating a variety of different land surface types, 
the final analysis consolidates all classes into four categories: — developed, forest, 
agriculture (including pasture and bare land) and water. A summary of these four classes 
by percent is listed below in Table 2. The table shows a progressive increase in developed 
land, which occurred largely at the expense of farmland. Unlike the corresponding values 
for State College, forested area in Chester Co. underwent a small but non-trivial decline. 
Much of this decrease in forested area pertains to the loss of trees in and around new 
developments, rather than the loss of land formally designated as forest. Instead, new 
housing in wooded areas in at least partly responsible for the decline in the forest 
category. 

Table 2: Percentage Classified for Each Category, Chester Co. 


Sfc. Type 

1987 

1988 

1991 

1993 

1996 

Agriculture 

71.8 

70.8 

69.3 

67.5 

66.0 

Forest 

16.8 

16.9 

15.1 

14.7 

14.7 

Developed 

10.9 

11.6 

15.1 

16.3 

19.0 

Water 

0.5 

0.6 

0.45 

0.38 

0.3 

Cloud, etc. 

0.01 

0.15 

0.0 

1.02 

0.0 
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Validation 


To verify these statistics, several approaches have been taken. The first was to compare 
our designated categories with those determined by aircraft and in situ methods, as 
published by the Chester Co. planning office in West Chester, PA (courtesy of Mr. Wayne 
Clapp), with whom we have maintained a cordial relationship. The published classification 
for 1987 was compared with our 1987 classification based on the Landsat image. We then 
determined a verification matrix from which we calculated user accuracy and producer 
accuracy (Congalton, 1991). These results are shown in Table 3. 

t 

I 

Table 3. User and Producer Accuracy: Chester Co. (1987 Survey) 


Sfc. Type 

User 

Producer 


Accuracy 

Accuracy 

veg/crops 

76 

96 

forest 

87 

75 

developed 

93 

54 

water 

100 

100 


User accuracy is calculated by dividing the number of correctly classified pixels in a land 
cover category by the total number of pixels classified in that particular category. This 
statistic represents the probability that a pixel classified using remotely sensed data 
actually represents that class. The producer’s accuracy is calculated by dividing the 
number of correctly classified pixels in a land cover category by the total sum of that 
category. This represents the probability that a reference pixel is correctly classified. 

For example, the producer of the 1987 classified image can claim only that 54% of the 
time an area that was actually developed was identified as such. However, a user of the 




classification map will find that if they were to visit a site that the map indicates is 
developed, 93% of the time it will have been developed. 

A second method, which is based on determining the so-called the KHAT statistic 
(Congalton, 1991), involves not only the correctly classified pixels but the off-diagonal 
elements in the matrix. KHAT takes into consideration confusion among classes. For the 
1987 image, the result is a KHAT of 81%, which implies that the classification process 
avoided about 4/5 of the errors that a completely random classification would generate. 
This result is in agreement with other remotely derived land cover classifications (Tin- 
Seong, 1995). 

Producer accuracy tends to be lower for developed areas than non-developed areas 
because of confusion between medium density buildings and bare soil; Tin-Seong (1995) 
found that this ambiguity accounted for about one-third of the error. In our 1987 image, 
developed areas were predominantly misclassified as either vegetated grassy field or 
forests. Spectral confusion of land cover types is aggravated by the spatial heterogeneity 
and mixed pixels typical of urban settings. Developed areas tend to include a considerable 
interpixel and intrapixel variation which is found less over forests or water bodies. Above 
all, a distinction must be made between spectral and societal classifications. Both may be 
correct within their own context. 

Additional significance testing for the Landsat classification when combined in 1 km 
squares ( ref. discussion below) in order to assess the confidence in determining the 
percentage of developed area. Classified Landsat TM images were subset into 1 km 
squares (roughly comparable to the resolution of the AVHRR), each containing 1 600 
pixels assigned to one of the five categories listed in Table 2. From this arrangement, a 
2x2 matrix was created for each of the 1 km pixels for two images at two different years; 
this is shown below in Table 4. 
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The diagonal elements represent pixels which have experienced no change in classification 
between the two times in question. Off-diagonal elements represent pixels which have 
changed from either undeveloped to developed or vice versa. A null hypothesis of no 
change over time is equivalent to a hypothesis of symmetry in the matrix. No development 
in time requires that the off-diagonals must be either zero or that X 12 = X 2J , which is to say 
that as many pixels changed from undeveloped to developed as from developed to 
undeveloped. 


Table 4. Matrix for obtaining p-value for development assessment 


Time 1 

Time 2 

Developed 

Not Developed 

Developed 

x„ 

X] 2 

Not Developed 

x 21 

x 22 


Following Eliasziw’s (1991) McNemar test for a 2x2 table, we are testing the null 
hypothesis of marginal symmetry by calculating the p-value: (p(% 2 >z 2 ), where 

_j_ {X n -Xnf 

(*I2 *21 ) 

The p-value is determined by comparing the z 2 value to a chi-square distribution density 
function with one degree of freedom, which yields the probability of obtaining a y2 greater 
than or equal to the calculated z 2 value. A low probability indicates strong evidence that 
the off-diagonals are not symmetric and there has been a significant change in land use 
over time. 

Qualitatively, at least, the general development trend demonstrated by Landsat TM 
classified images reflect the changes which have been reported for Chester Co. by the 
Chester Co. Planning Commission. These patterns of change, shown by the classified TM 
images for 1987 and 1996 (Figures la and b) reflect continued suburban expansion 
outward from Philadelphia (denoted by the letter ‘a’ in Figures la and b) into eastern 
Chester Co., as well as the development of ‘bedroom communities’ out of Delaware 
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June 10. 1987 
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Chester County, PA and Surrounding Area 


4460000 


4450000 


4440000 


4430000 


4420000 


4410000 


4400000 


420000 430000 440000 450000 460000 470000 480000 490000 



420000 430000 440000 450000 460000 470000 480000 490000 


4460000 


4450000 


4440000 


4430000 


4420000 


44 10000 


4400000 



Legend 

bare soil 1 

bare soil 2 

partially vegetated 

scrubby vegetation 

vegetated / grassy field 

forest 

water 

developed 

cloud cover 


Scale 


20 


3 Kilometers 
0 


10 


0 


Miles 



derived from Landsat TM imagery 
georeferenced to UTM map coordinates, Zone 18 


Figure la. Classified Landsat TM image for June 10, 1987 over Chester Co. Letters denote 
locations referred to in the text. The outline of Chester Co. is shown by a thin solid border in the 
background. Categories referred to in the text (Tables 2 and 3) correspond to those in the legend 
except that the agricultural classification includes the five bare soil and vegetation classes. 
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June 2, 1396 

Land Cover Classification 
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Figure lb. Same as Figure la except for 2 June 1996. 
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(location ‘e’ in Figures la and b). At a first glance, we are likewise finding the most 
development between 1987 and 1996 in the West Chester borough and the surrounding 
townships along route US 202 in the east and along US 30 which runs through the middle 
of Chester Valley (location ‘s’ in Figures la and b). 

Present development, as revealed by Landsat, presage the Planning Commission’s 
projected development trends for the year 2000, in which the eastern part of the county 
along US 30 becomes densely development and the previously distinct centers of 
development become more of a continuous mass stretching through the eastern half of the 
county. Although major roads serve as conduits for rapid and dense urbanization, many 
tiny, isolated pockets of development are already starting to form in the highly rural 
southwestern part of the county (location ‘t’ in Figures la and b). These areas are 
projected to increase in size and number. 

The decrease in forest from 1987to 1993 reflects this development pattern, as most of the 
land use change occurred along the US 30-202 axes. In the classification, the highly 
wooded residential areas of this region often appear as grids of forest and developed 
classes. With continual expansion of the urban fringe out of Philadelphia, the developed 
class expanded significantly in 1996, resulting in a loss of forest class. No loss in forest 
was detected between 1993 and the present. 

We have also attempted to compare the 1996 classification with in situ observations made 
in the field. In September, we conducted a field survey of Chester Co. during which time 
we were able to locate (using a Trimble™ GPS meter) about 45 pixels in the Landsat 
image and to classify the locations visually. These statistics are currently being processed. 

A VHRR images 

AVHRR imagery is now being resampled and superimposed on the 1 km squares used to 
create the z 2 statistics referred to in Table 4. Preliminary results tend to support those of 
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Owen and Carlson (1997) which show that (1) neighborhood scale changes, referred to 
above, can be detected by A VHRR when changes in developed area is greater than 1 0%; 
(2) vegetation, notably tree cover, affects Mo more than individual precipitation events or 
protracted dry or wet periods; Mo appears to be an intrinsic descriptor of land surface; (3) 
Mo responded very little to changes in development for a given type of development; (4) 
fractional vegetation cover and scaled surface radiant temperature changed significantly in 
response to urbanization; (5) evapotranspiration decreases at approximately the same rate 
as the increase in developed surface area. 

At present, Chester Co. images are begin transformed to values of Mo, Fr and 
evapotranspiration and their values at fixed UTM locations will be charted at selected 
locations over the past 10 years. As in the accompanying paper by Carlson and Owen 
(1996), pixel migrations will be followed within the context of the triangle. We plan to 
obtain additional land use and other socio-economic data from the Chester Co. planning 
office. Scatterplots and analyses of Mo, Fr and evapotranspiration should be finished by 
early 1997. 

An example of the land use classification is shown in Figures la and lb. 

Costa Rica 

An allied project involving deforestation in Costa Rica. Deforestation has been 
considerable in this country during the past 30 years. Initially, our approach will follow 
procedures developed for Chester Co. Unlike Chester Co and State College, land use 
changes have come largely at the expense of forest land, rather than agriculture. A new 
and vital element in our Costa Rican program is a collaboration with Professor Arturo 
Sanchez (currently at the University of Costa Rica), who has acquired and georeferenced a 
number of Landsat TM images. Professor Sanchez brings to the project a considerable 
knowledge of the forest resources, ecology, hydrology and demographics of Costa Rica, 
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as well as his expertise in image processing and analysis. He has already prepared several 
extensive data bases for Costa Rica. 


To this data set provided by Dr. Sanchez, we have added 5 more Landsat images and six 
AVHRR images over a 10 year period, 1986 - 1996. AVHRR scenes will be 
georeferenced to the Landsat images, which have been rectified by Dr. Sanchez. A 
summary of these images is presented in Table 5. A few more images for 1995 and 1996 
may be collected before the end of the project. 

Table 5: List of Scenes Obtained for Costa Rica; 
month and year). (AVHRR and Landsat TM) 


Landsat-TM 

Landsat-TM 

AVHRR 

3/75 

09/88 

12/85 

01/76** 

02/89 

12/90 

03/78 

12/90** 

01/92 

12/78 

03/91 

12/92 

01/79 

03/91 

12/93 

01/84 

12/91 

| 01/95 

02/86** 

03/92 

| 03/96 

02/86 

09/92 


08/86 

03/93 

t 

01/87 

03/93 


01/87 

04/93** 

; 

11/87 

03/94 


12/87 

03/94 

< 


03/96** 



** Indicate usable images for San Jose area 
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The study centers on the capital area, San Jose’ and its surroundings, where the major 
effects on the forest and the countryside have occurred. Not all images are usable because 
of extensive cloud contamination, although all dates correspond to the dry season, which 
extends from December through March. Starred entries in Table 5 refer to Landsat images 
which are sufficiently cloud free for the analysis of the deforestation and urbanization in 
the target area. 

Georeferencing of TM images is being performed at the University of Costa Rica by Dr. 
Sanchez and later this winter and at Penn State for AVHRR. The PI will make one trip to 
Costa Rica in January and Dr. Sanchez is expected to pay a working visit to Penn State in 
December. Like our Chester Co. project, analyses will be made within the context of the 
triangle and the results related to socio-economic and ecological changes. In regard to the 
latter, we will attempt to use the changes in Mo, Fr and evapotranspiration to describe 
how the deforestation exerts pressures on the boundaries of ‘life zone’ (regions of uniform 
ecology) in Costa Rica. This part of the project will be outlined in more detail in a later 
report. 

3. Operational Goals 

Not only does scaling allow us to compare pixel values from differing images but the 
effect of scaling greatly minimizes errors due to atmospheric correction, instrument drift 
and viewing angle and therefore reduces the importance of obtaining absolute accuracy in 
deriving the surface parameters. Indeed, preliminary results using radiative transfer model 
simulations (Carlson and Ripley, 1997) shows that scaling ND VI and surface radiant 
temperature yields almost identical values of fractional vegetation cover and surface 
moisture availability as would be obtained from radiances corrected with MODTRAN, a 
recent variant of LOWTRAN (Kneizys et al., 1988). 
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This is illustrated in Figure 2, which shows scaled NDVI ( N* ) based on the uncorrected 
(at sensor; N* 2 (a)) and corrected (at surface; N* 2 (c)) radiances. The figure not only 
confirms two earlier and independent validations of a square root law between scaled 
NDVI and fractional vegetation cover (Choudhury et al., 1994’ Gillies and Carlson, 
1995), but suggests that the relationship holds equally well for uncorrected as for 
corrected radiances. 



Figure 2. Simulations with a simple radiative transfer model for 
scaled NDVI squared (N* 2 ) versus fractional vegetation cover (Fr) 
based on the uncorrected (at sensor; N* 2 (a)) and corrected (at 
surface; N* 2 (c)) radiances. The solid sloping line without symbols 
represents a 1 : 1 correspondence between N* 2 and Fr. 


Scaling of temperature and NDVI not only facilitates comparison between images and 
scatterplots, but it greatly reduces measurement error, eliminating the need for a detailed 
knowledge of plant, soil or atmospheric conditions. We can also show that scaling 
temperature greatly reduces the atmospheric correction to the thermal radiances, just as 
scaling NDVI virtually eliminates the need to correct N* for atmospheric attenuation 
(Figure 2). 


To the extent that the uncorrected and corrected NDVI values are linearly related, the 
atmospheric correction is eliminated from N* in scaling. Indeed, we find this to be the case 
(Carlson and Ripley, 1997). The result of scaling is that virtually identical values of Mo 
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and Fr are obtained for corrected and uncorrected values of NDV1. Errors in Mo 
produced by not scaling the surface radiant temperature occur because the uncorrected 
temperatures are not linearly related to the corrected surface temperature. In split window 
models, such as that of Price (1984), the atmospheric corrections depend on differences in 
apparent temperatures between channel 4 and channel 5 of AVHRR, and so are not 
necessarily linearly related to the temperature of either channel. 

Consider first a scatterplot for one of the Chester images (Figure 3). Using the Price split 
window method we first corrected every surface temperature and then found a linear 
regression between the temperature corrected with MODTRAN for channel 4 of AVHRR 
and the apparent (uncorrected) temperature, as measured by the satellite. This established 
a unique correction for each value of surface temperature although the slope and intercept 
of the regression are based on the split window results. RMS differences between this 
straight line fit to the split window temperatures and the temperatures given by the split 
window algorithm for each pixel provide an estimate of the ‘typical error’ one might 
encounter by collapsing the split window algorithm to a simple linear regression. (The 
word ‘error’ here is simply a measure of the scatter from a straight line, as there is some 
question as to whether this scatter represents noise or real variations.) The RMS error in 
the (non-dimensional) normalized temperature (T*) for the test case (Figure 3) was ±0.08. 
This value constitutes a relatively small fraction of the full range in T*, which varies from 
zero to one, respectively, between air temperature and the radiometric temperature of dry, 
bare soil. 

To obtain isopleths of Mo, a polynomial relationship between Mo, N*, and T* is 
determined from making a range of simulations with our 

soil/vegetation/atmosphere/transfer (SVAT) model (Gillies and Carlson, 1995). From this 
output, Mo is expressed as a function of the two scaled variables; isopleths of Mo are 
shown in Figure 3. 

From this we also compute the following expression: 
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AM 0 (N*, T*)*cM 0 /?T* [AT] 


where the partial derivative is determined by the functional form of Mo lines obtained 
from the simulations (the sloping lines in Figure 3) and the term in brackets on the right 
hand side is the RMS error referred to above. Note that this is really an estimate of the 
expected error in Mo based on a single case and may not precisely represent the error one 
would obtain in not using the split window on any other case. Further analyses of the 
remaining 13 AVHRR images for the Chester Co. area is being undertaken. 

The implications of Figures 3 and 4 are similar to those for Figure 2 in the case of N*. The 
residual error in not making atmospheric corrections to surface radiant temperature is 
generally less than 0.2 in Mo (20% of the full range) except for the ranges of large Mo 
(above 0.5 - 0.8) and large Fr (above 80%). Inasmuch as most of the points in this 
scatterplot lie outside this range (except for densely vegetated areas) in this test case 
(Figure 3), the percentage of pixels suffering a significant degradation in Mo due to 
neglect of an atmospheric correction is quite small. Two important implications can be 
drawn from this result: First, scaling is a powerful tool which may permit the use of 
uncorrected surface radiant temperatures, at least in an operational mode, to determine 
Mo and Fr from satellite data. This shortcut would significantly reduce the work required 
for data reduction. Second, use of a split window algorithm may not yield significant 
improvement over a linear least squares relationship between corrected and measured 
surface radiant temperature. 

It is our intention to make the triangle method useful to the scientific community, which is 
to say for use by other scientists. The methodology is currently being prepared for display 
on our web site. As a first step we are creating a web site (http://www.essc.psu.edu/~tnc), 
from which one can access instructions on how to apply the triangle method and how to 
access and execute the algorithms from which Mo/Fr isopleths and scatterplots can be 
generated. 
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Figure 3. Scatterplot for the NOAA/AVHH image over Chester Co., 25 
July, 1985. Normalized surface radiant temperature (T*) versus scaled 
NDVI (N*). The sloping vertical lines are isopleths of moisture availability 
(Mo) 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 


Normalised Temperature 

Figure 4. Normalized surface radiant temperature (T*) versus scaled NDVI 
(N*). The shaded area corresponds to the domain of the shaded region (the 
triangle) in figure 3 and the dark shaded area to the part of the triangle 
where the RMS ‘error’ in Mo (the RMS difference between Mo obtained 
from surface radiant temperatures corrected with a split window algorithm 
and those from uncorrected temperatures) exceeds 20% of the value of Mo. 
(Based on the NOAA/AVHH image over Chester Co., 25 July, 1985.) 
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4. Summary and Future Work 


We intend to extend our State College analyses in new directions, tying the land use 
changes within the context of the triangle to socio-economic changes in the region and to 
possible local climate changes and changes in surface hydrology produced by urbanization. 
Ideally, we would like to produce a description of land use changes and urban 
development for Chester Co. in terms of changes in the climate parameters, Mo and Fr, 
and in the land use index developed by Owen and Carlson (1996). For Costa Rica, we will 
relate the changes in land surface parameters, Mo and Fr, to the so-called ‘life zone’ which 
represent areas of more or less uniform ecology. Although most of the rampant 
deforestation in Costa Rica has been checked during the past 5 years, our image library 
extends back to the 1970s when deforestation was among the highest in the world. 

Preliminary results for Chester Co. show that urbanization can be detected throughout the 
county during the period 1987-1996, with an increase in percent of area developed from 
about 1 1 in 1987 to about 19% in 1996. Urbanization is exploding throughout the 
countryside and expanding rapidly from roads and highways. 

Chester Co. analyses should be completed early in 1997; all Costa Rican scenes will be 
completed and scatterplots prepared by mid 1997. Although the Chester Co. work should 
be completed by the end of summer 1997, we are hopeful that this research will continue 
in some form beyond 1997. In anticipation that the Chester Co. results will prove to be 
very interesting, we have requested continued satellite coverage of both Chester Co. and 
Costa Rica from the ASTER project, beginning in 1998. Coverage requests, recently 
submitted to the ASTER project at NASA, are modest; they amount to two or three 
scenes every summer for Chester Co. and about the same number of scenes during the dry 
season (Dec - Feb) for Costa Rica. ASTER’s 30 m surface resolution will allow it to 
replace Landsat for image classification, but its several thermal IR channels will also 
permit the use of split window corrections schemes for surface radiant temperature, if 
indeed these corrections are required. This is currently possible only with AVHRR. 
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Depending on the morning overpass time, ASTER may also replace AVHRR, whose 
overpass times are close to that of maximum surface temperature. 

It is our hope that the triangle method can be simplified and streamlined such that it will 
lend itself to a routine applications. In so doing, it may be possible to use the fractional 
vegetation cover and moisture availability not only as key land use parameters in climate 
models but as unique indices of land use change which provide added descriptors of 
urbanization and deforestation beyond the conventional parameters already being 
employed by county planning boards. 

Finally, we will be presenting a poster at an international meeting on remote sensing in 
France, (Physical Measurements & Signatures in Remote Sensing), to be held in 
Courchevel between 7 and 1 1 April, 1997. The subject of our poster is urbanization and 
deforestation in Chester Co. and Costa Rica, as viewed using the triangle method. 
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PURPOSE: 

Urbanization — the creation of industrial, commercial and residential spaces from agricultural or forested 
land — has attracted the attention of climatologists since the early part of this century when the urban heat 
island came to be well documented. Although much has been written about the impact of urbanization on 
surface energy budgets and on surface air and radiant ground temperature, difficulties remain in translating 
this information into practical applications. We see this problems in three parts: (1) the need to translate 
changes in urban land use into fundamental state parameters (other than temperature) which could be used, 
given a particular land use change, to both describe and predict the resulting changes in urban 
microclimate; (2) urban planners have not yet felt the need to make use of potential climate descriptors, and 
(3) remote sensing is not yet appreciated as a tool in monitoring neighborhood-scale land use changes over 
periods of years or decades. The purpose of this paper is to indirectly address these three issues by 
illustrating a method for analyzing the climatic impact of urbanization using satellite measurements. 


1. OVERVIEW 

Unrestrained growth of cities, depletion of natural 
resources, the disappearance of forests along with 
Midlife habitats are among the main concerns of 
humankind as the 20th century draws to a close. 
Much of recent urban expansion in the United 
States has taken place at the expense of agricultural 
land, but significant loss of forested areas has also 
accompanied the growth of outlying urban centers. 
Such a transformation in the land surface 
profoundly affects the local microclimate. 

However, not all such changes in microclimate are 
necessarily deleterious to humans, nor is 
urbanization always accompanied by an increase in 
temperature and a decrease in vegetation amount 
and soil moisture; the opposite may be true for 
desert environments. 

Our ability to monitor the climatic effects of either 
the beneficial or deleterious effects of deforestation 
and urbanization is hampered by the scale of the 
encroachment, by the inaccessibility of many 
affected areas, by constraints imposed by regional 
political sensibilities, and by insufficient interest or 
insufficient funds in the scientific community. 
Nevertheless, urbanization is recognizable on 


satellite images by changes in surface radiant 
temperature and vegetation index, from which 
estimates of soil moisture, surface runoff and 
vegetation fraction can be made. Most of the 
attention in satellite remote monitoring of 
urbanization has been focused on the so-called 
'urban heat island'. It has become increasingly 
clear, however, that urbanization no longer involves 
a simple expansion or development of a 
recognizable population center but, rather, the 
formation and expansion of many separate 
commercial and residential centers about a 
relatively stable inner region. 

Relating changes in local microclimate due to 
urbanization to the impact it has on human activity, 
such as the added or reduced energy consumption 
in buildings, changes in the length of growing 
season, etc., still remains an elusive goal. More 
important, however, than deriving a single set of 
parameters' is a description of the relative and 
ongoing changes in microclimate that occur over 
years or decades and the ability to predict future 
changes given a specific growth scenario expressed 
in terms of changes in derived land surface 
parameters. 



With very few exceptions, land use classification 
data, as represented in Geographical Information 
Systems (G IS), arc expressed in terms of non- 
physical descriptors, such as commercial, 
industrial, forest, or agricultural. For example, 
such information is the basis of detailed maps of 
deforestation (e.g., Sader and Joyce, 1988). 
Furthermore, many investigations of deforestation 
and land use change, based on remote surface 
measurements of land use and vegetation cover, 
tend to focus on a general description of the 
vegetation, i.e. , on estimating biomass or green 
foliage (Tucker et al, 1984; Asrar et al, 1984; 
Goward et a /. , 1991), or as a qualitative tool to 
study specific ecosystems (Hall et al ., 1991). Sader 
et al. (1990), in an overview of current efforts to 
apply remote sensing techniques to monitor tropical 
forests, merely states that current satellite 
technology is adequate for this purpose. Thermal 
infrared data, however, is usually not incorporated 
in this analysis, but is relegated to a secondary role, 
such as to assess the degree of cloud 
contamination. 

These semi-quantitative methods are, of course, 
highly useful for a variety of purposes (sociological 
and ecological) but lack the physical basis which 
permits the assessment of changes in microclimate. 
Moreover, attempts to monitor deforestation over 
extended periods have proved elusive, and are often 
highly qualitative (Sader et al. , 1990), although 
Tucker et al (1985) express a guarded affirmation 
that satellite measurements of vegetation indices 
reflect trends in agricultural practices. 

Remote measurements of thermal infrared 
temperatures have been used much more widely in 
studies of urbanization than for deforestation. In 
most cases, the focus has been on the temperature 
perturbation known as the urban heat island 
(Matson, et al, 1979; Goward et al, 1985; Roth et 
al y 1989). Recent investigations, however, have 
begun to incorporate vegetation index 
measurements, along with surface radiant 
temperature differences between urban centers and 
their surroundings, as a measure of urbanization 
(Gallo et al, 1 993). While such papers tend to 
show a relationship between the urban-rural 


differences in vegetation and surface radiant 
temperature, they provide no information on the 
fundamental parameters from which predictions of 
future surface climate changes can be made. 

There is, however, no doubt that a possibility exists 
to estimate, with the aid of satellites, surface 
variables that reflect the impact of deforestation 
and urbanization from the microscale 
(microclimate) to the macroscale (regional) climate. 
For example, there are clear relationships between 
variables such as water stress and vegetation 
biomass as well as a number of radiance-derived 
indices (e.g., Nelson e/a/., 1987; Szekielda, 1988; 
Pierce et al. , 1 990; Spanner et al , 1 990; Cihlar et 
al. , 1991; and Hall et al 1991). Nemani and 
Running (1989) suggest a relationship between 
normalized difference vegetation index (NDVI) and 
a canopy resistance for deriving a measure of water 
stress in forests; while Thunnissen and 
Nieuwenhuis (1990) and Jackson et al (1977) 
discuss an empirical land surface parameter 'B' (a 
constant), from which one obtains daily surface 
evapotranspiration. Carlson et al (1981) also 
derived two bulk land surface parameters (the 
thermal inertia and moisture availability) 
representative of the potential for 
evapotranspiration at the surface. Many of these 
relationships, however, are based on bulk 
parameters, i.c. y those which represent the 
combined impact of soil and vegetation. 

2. PHYSICAL BASIS 

Current state of the art boundary layer schemes 
(e.g., Taconet et al, 1986; Sellers et al , 1986, 
Abramopoulos et al y 1989; Dickinson et al , 1991) 
demonstrate the importance of calculating separate 
energy budgets for both soil and vegetation 
surfaces (expressed as proportions of soil and 
fractional vegetation cover). This would make 
surface moisture availability and fractional 
vegetation cover fundamental parameters in 
determining the surface energy budget (Avissar, 

1 992; Koster and Suarez, 1 992; Carlson et al , 
1994). 



Carlson el al (1994) define the surface moisture 
availability (M 0 ) as the ratio of extractable surface 

soil water content to field capacity and the 
fractional vegetation cover (Fr) as the fraction of 
the surface visibly covered by leaves. As 
demonstrated by Collins and Avissar (1994) and 
Friedl and Davis (1994), under clear sky conditions 
and during day-time hours, the surface 
evapotranspiration is largely a function of surface 
moisture and vegetation cover. Thus, both 
parameters have precise physical meaning and are 
the important factors in driving the surface energy 
budget. As we shall illustrate, they constitute a 
Land Cover Index (LCI) which can be used to 
describe urbanization and deforestation, and serve 
as initial boundary conditions in the land surface 
components of atmospheric models (Smith el al 
1994). The LCI should respond to normal seasonal 
evolution of the vegetation canopy or, more 
importantly, to direct human impact on the land 
surface. 
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where ^ is the specific humidity of the air, q s (T s ) 

the saturation specific humidity at the surface 
radiant temperature T ir , Le the latent heat of 

vaporization and p the air density. 



2.1 Some Definitions 

Isopleths of the surface moisture availability in the 
schematic Figure 1 are derived using the triangle 
method, so called because the distribution of pixels 
in surface temperature/vegetation space tends to 
form a triangular envelope (Carlson el ai y 1994; 
Gillies and Carlson, 1995; Gillies el al,. 1996). 
This method uses normalized difference vegetation 
index (NDVI) and surface radiant temperature (T ir ) 

to obtain surface soil water content and surface 
energy fluxes. NDVI, which is equated with 
fraction vegetation cover, is defined as 

NDVI — a ~) . (!) 

(a . + a . ) 

where a represents reflectance and subscripts nir 
and vis refer to wavelengths above and below 0.7 
microns. Surface moisture availability is defined in 
terms of the surface evaporation over bare soil, 
which is: 


Figure 1. Schematic illustration of the triangle 
showing the variation of vegetation cover, surface soil 
moisture availability and surface radiant temperature. 
The darkened right border denotes the ‘warm edge* 
and the arrow indicates a pixel migration in response 
to urbanization. 

In order to compare scatterplots from different 
images, NDVI is linearly scaled between limits for 
bare soil (NDVI 0 ) and 100% vegetation cover 

(NDVI S ) (Choudhury el al. 1994; Gillies and 

Carlson, 1995). Scaling of NDVI is necessary 
when images from several days are to be compared 
because the absolute values of NDVI tend to vary 
temporally in a non -systematic manner (Price, 
1987). Therefore, the scaled value is defined as 

(NDVI - NDVI e ) 

(NDVI,-NSV1 0 ) U 

Temperature is also scaled between the warmest 
(M o =0; NDV7 0 ) (T^) and coldest (M 0 =l; NDVI S ) 

(T») comers of the triangle, as indicated in Figure 



2; these are, respectively dry', bare soil and wet soil 
at 100% vegetation cover according. A scaled 

temperature ( T) is similarly defined as 
(T(Mo = Q\NDVI )-T ir ) 

*j> o 

” ( T(Mo = 0 \NDVI o )-T(Mo = 1; NDVI S )) 

(4) 

where T varies from zero to one as does M 0} whose 
isopleths extend from vertex to base of the triangle 
These limits are determined first by inspection of 
the unsealed scatterplot supplemented by 
simulations with a 

soil/vegetation/atmosphere/transfer (SVAT) model, 
from which isopleths of M 0 and the surface energy 
fluxes are obtained. As the cold anchor point 
(M 0 =l; NDVIj) is consistently close to air 
temperature, it is customarily assigned the ambient 
air temperature (T»). 

An assumption in the triangle method is that M 0 
varies monotonically at constant NDVI from zero 
on the right side (the ‘warm edge’) to 100% of its 
maximum on the left side and NDVI varies 
monotonically with Fr. We find that the empirical 
relationship, Fr «N* 2 , applies to a number of case 
studies (e.g. t Gillies et al. 9 1996). 

The importance of Fr in the surface energy balance 
can be seen through the following equation for the 
total evapotranspiration (LEu*), which depends on 
the evaporation from bare soil (Equation 2) and 
transpiration from the leaves (LEf) 

LE lol =LE 0 (\~Fr) + LE f {Fr) (5) 

Gillies et al (1996) argue that LEf is usually close 
to the potential (maximum) value for a particular 
set of meteorological conditions unless the plant is 
undergoing severe water stress. As the temperature 
scaling (Equation 4) tends to remove fluctuations 
produced by changes in ambient temperature, both 
spatial and temporal variations in the surface fluxes 
and surface radiant temperature over a mix of 
vegetation and bare soil are strongly modulated by 
the fractional vegetation cover, Fr, and by M 0 via 


Equation 2. However, the sensitivity of the latent 
heat flux to M 0 and Fr occurs largely below about 
0.4 (or 40%) (Figure 3). 



Surface Radiant Temperature { °C) 


Figure 2. Scatterplot of surface radiant temperature 
versus NDVI over Center County, PA, for an image 
during June of 1994, with the warm edge labeled. 
Anchor values for scaling the triangle are indicated: 

T i( observed air temperature, T^, the temperature of 
bare soil and zero extractable soil water (Mo=0), bare 
soil (NDVI 0 ) and 100% vegetation cover (NDVI,). 

3. ESTIMATING Mo AND Fr: AN 
ILLUSTRATION 

3.1 Methodology 

Gillies and Carlson (1995) demonstrate that 
urbanization can be viewed in the context of Figure 
1, in which a surface location undergoing 
urbanization over time will usually experience an 
increase in T^ resulting from a reduction in both Fr 
and M 0 . The arrow in Figure 1 illustrates this 
hypothesis. Relating land cover statistics with 
values of Fr and M<> requires the coupling of sub- 
pixel urban land cover data to pixel-wide values of 
Ty and NDVI. The former was attained through the 
Landsat TM sensor and the latter through the 
National Oceanic and Atmospheric Administration 
(NOAA) Advanced Very High Resolution 
Radiometer (AVHRR). 





Figure 3. Sensitivity of latent heat flux LE tot (Win* 2 ) to 
Mo and Fr, as obtained from simulations with the 
SVAT model. 

Figure 4 shows the location of the research study 
area, which encompasses a 1000 km 2 region of 
Centre County, Pennsylvania, USA and includes 
the borough of State College. This study area was 
selected due to its limited size, tight land-use 
zoning, isolation from other metropolitan areas and 
convenient access for field verification of the 
results. 

With a growth rate of 9.8% during the 1980s, 
Centre County was a site where substantial urban 
development had taken place over the period in 
question (Centre County Planning Office, 1992). In 
light of this, eight AVHRR (in Level lb format) 
scenes were obtained for the study area from the 
U.S. Geological Survey's EROS Data Center- two 
single-date. Local Area Coverage (LAC) scenes 
each for the summers of 1985, 1988, 1991 and 
1994. These scenes were selected to minimize both 
cloudiness and large viewing angles from nadir of 
the satellite; viewing angles were restricted to less 
than 40° to minimize unacceptable errors in 
reflectance values in lieu of performing a rigorous 
correction for angular dependency of reflectance 
(Paltridge and Mitchell, 1990). Given these criteria, 
coupled with the requirement that images be 
selected during the period of stable vegetation 
growth (from mid-June through mid-August) 

(Fisher, 1994), a limited number of scenes were 
available for selection. 



Figure 4. Location of the target area: Southern Centre 
CO, PA. 

It was not possible, however, to completely avoid 
cloud and haze contaminated pixels in some scenes. 
Fortunately, such pixels were identified and 
excluded from the data analysis by using a 

modified cloud screening technique of Carlson*?/ 

al. (1994). Once selected, the raw data for each of 
the eight AVHRR scenes was rectified and 
radiometrically corrected in order to assure a 
precise spatial and temporal inter-scene comparison 
(Emery et al 1989). Two Landsat TM scenes 
were used in this study— one from 1986 and 
another from 1993. Four pairs of AVHRR images 
were analyzed, two for each of the summers of 
1985, 1988, 1991, and 1994. Image pairs were 
combined into a single scatterplot for each of the 
four years. All AVHRR images pertain to relatively 
undisturbed days not immediately following 
significant rainfall. 

A subset of each of the Landsat TM images (1986 
and 1993) was taken to encompass the study area 
(Figure 4) and subsequently rectified (resampled to 
25 meters resolution with a nearest neighbor 




algorithm) to UTM (Universal Transverse 
Mercator) coordinates. This method of resampling 
was also applied to the AVHRR data as 
maintaining the different and discrete levels of 
temperature was of primary importance while in 
the case of the Landsat TM data set nearest 
neighbor resampling is considered sufficient before 
any form of classification is undertaken. 

Radiometric correction was performed on each of 
the AVHRR scenes by first applying calibration 
coefficients provided in the ephemeris data, 
followed by additional algorithms. In computing 
surface reflectance (/.<?., correcting for the effects 
of atmospheric scattering etc.) value for channels 1 
and 2, a bulk haze correction was applied, based 
upon an offset that corresponded to the minimum 
reflectance value for each channel (Richards, 

1993). The at-sensor radiances for channels 4 and 
5 were converted to at-surface-radiant 
temperatures. 

To facilitate an accurate temporal comparison 
between NDVI values from all AVHRR scenes, a 
sensor degradation correction was applied to all 
NDVI values, based on the number of months 
between satellite launch and scene overpass (Che 
and Price, 1992). This corrected set of NDVI 
values was the final product used in the data 
analysis. 

Land cover statistics were derived for both Landsat 
TM scenes (1986 and 1993). These scenes were 
initially subjected to a terrain normalization 
procedure that minimized solar illumination 
differences caused by topography (Hodgson and 
Shelley, 1994). Then, a maximum likelihood 
supervised level I classification of developed areas 
was performed according to Anderson el ai 
(1976). Five classes were distinguished— water, 
forest, agriculture, developed and unclassified and 
ground-truthed using the following— a 1990 EPA 
land cover map, Centre County aerial photographs 
and expert knowledge from Centre County's 
Planning Office (Pennick 1995, private 
communication). The confidence of each pixel 
classification was evaluated using a spectral 
distribution of pixel variances from the mean of 


each class. Those pixels that exceeded a confidence 
interval threshold of 95% (85% for urban) using a 
Chi-square statistic were reclassified as unknown. 
From the five known classes, land cover statistics, 
that included the percentage of urban development 
for each AVHRR pixel, were compiled from a 
possible 1600 Landsat TM pixels corresponding to 
each AVHRR pixel. 

3.2 Results 

Land cover was classified into three broad classes: 
forested, agricultural and developed. AVHRR 
pixel subsets for each of these classes were based 

upon the classified Landsat TM images. T and Fr 
attributes for each class were distinctly separate. 
The forested pixel average, for example, was 

centered near the point T = 0.17, Fr = 0.8, while 
agricultural and urban pixel averages were centered 

at higher values of T and lower values of Fr. The 
separation of urban areas from other land cover 

types in terms of T and Fr reflects the sensitivity 
of vegetation and surface radiant temperature to 
urban land use. High spatial resolution thermal 
infrared measurements by Nichol (1996) suggest 
that such distinctions, particularly between urban 
and vegetated pixels, yield similar separation in the 
surface radiant temperature response at the 
resolution of a city block. Although no single pixel 
in the study area was completely urbanized, by 
implication, fully urbanized surfaces would be 

found where T > 0.9 and Fr = 0 beyond the lower 
right comer of Figure 5. 

From the distribution of land cover types in 

Fr/ T space, a land cover index (LCI) is determined. 
This allows the temporal process of urbanization to 
be quantified in terms of a location's starting point 
in the figure (representative of the mean influence 
of surrounding land cover) and the magnitude of 
change (representative of the urbanization process, 
assuming that surrounding land cover remains 
relatively unchanged during the time period under 
consideration). The spatial distribution of LCI 
values in the figure are determined by tangents to 
the one and two standard deviation ellipses normal 



to an overall linear relationship between Fr and T 
for all eight AVHRR images. The thin curved lines 
are isopleths of M c . 

The LCI and M 0 isopleths show the temporal 
relationship between the land cover parameters and 
urban land cover for pixels in the study area that 
underwent "robust development" (defined for pixels 
that underwent a minimum of 10% increase in 
development from 1986 to 1993 as defined by the 
two Landsat land cover classifications). Altogether 
5 areas exhibited more than 10% or more 
urbanization, as defined in the Landsat 
classification; (Two of these are shown in Figure 

5). 

For comparison, a composite trajectory of 
unchanged urban pixels is shown; this was derived 
from the mean of an ensemble of pixels that were at 
least 25% developed in 1986 and underwent no 
further increase in development by 1993. By 
definition, this control composite pixel should 
remain stationary in the absence of inter-scene 
variability. For reasons to be discussed in the next 
section, the 1988 points were excluded from Figure 
5. 

An interesting contrast, to be discussed in the next 
section, is between the Park Forest (10% increase 
in urbanization) and the Glenview (25% increase) 
neighborhoods. The difference in the 1985 Fr and 
M 0 values between the Park Forest and Glenview' 
pixels in Figure 5 is related to the dominant 
surrounding land cover class for each pixel. Park 
Forest is predominantly deciduous forest, w'hile the 
Glenview pixel consisted of formerly agricultural 
land with small homes, short or sparse vegetation 
and partially bulldozed fields. A comparison of the 
1986 and 1993 Landsat classifications shows that 
the make-up of surrounding land cover classes in 
these and all five analyzed pixels remained 
unchanged. Thus, urban development was the 
primary cause of the migration of the pixels in the 

Fr/ T space. 



Figure 5. Temporal migrations of pixels that 
experienced at least 10% development between 1986 
and 1993 (excluding 1988 data). The control trajectory 
refers to the mean values for a sample of pixels with 
over 25% development in 1986 that underwent no 
additional urbanization between 1986 and 1993. 

4. DISCUSSION AND CONCLUSIONS. . 

The Glenview and Park Forest trajectories, as well 
as those of the other three pixels (not shown) 
exceeding 10% urbanization, experienced similar 

decreases in Fr in conjunction with increases in 7\ 
and similar decreases in the LCI. Surprisingly, 
however,. M 0 exhibited little temporal variation. 
Although it may first appear that M 0 is low' at both 
sites, small differences in M 0 may be highly 
important in classifying a particular site, as can be 
seen in Figure 3. 

Another unexpected result was that the 1988 data 
(not shown), which corresponds to on of the 
warmest and driest summers on record at State 
College, exhibited little change in M 0 from the 
other years. In contrast, Fr decreased notably from 
1985 to 1988 to values near or below those of 
1994. at all five sites. Inclusion of the 1988 data 
therefore gives the trajectories a more zigzag 
appearance, but with an unmistakable decline in Fr 

and increase in T at nearly constant M 0 from 1986 
to 1994. 





Figure 6a. Photograph of Park Forest development 


Differences in M 0 between these Park Forest and 
Glenview sites are probably related to the tree 
cover and to the differences in the amount of bare 
soil visible to the radiometer. The striking contrast 
between Figures 6a and 6b suggests that Mo 
constitutes an intrinsic surface parameter 
independent of Fr. As Nichol (1996) has shown, 
bare soil can exhibit a higher radiant temperature 
than buildings, so that pixels closest to the warm 
edge may be more characteristic of undeveloped 
urban surfaces or unirrigated farmland than of 
center cities. 

Finally, a significant finding is that AVHRR data 
can be used to classify' the growth of individual 
housing developments. The triangle method 
provides a framework for examining this growth 
because it is relatively easy to apply and it is 
subject to a minimum of fluctuation due to sensor 
drift or changing atmospheric conditions. In 
portraying urbanization so graphically and in terms 
of parameters that govern the surface energy 
budget, the method may heighten awareness of the 
urbanization process . 



Figure 6b. Photograph of Glenview development 


The importance of trees in mitigating the climatic 
effects of urbanization - an idea which has begun 
to take hold in Atlanta, GA (Courtney, 1996; 
private communication) -- can not be understated. 
Values of LE lot calculated with the SVAT model 
indicate that Park Forest experienced 
approximately a 15% decline in evapotranspiration 
and Glenview a 30% decline over the nine-year 
period with values at Park Forest being much 
greater than those at Glenview. 

Dow and DcWalle (1995) show that urbanization 
produces a decline in evapotranspiration and a 
corresponding increase in runoff. If so, the planting 
of trees in urban areas would help to retain the 
precipitation in the ground, while not constituting 
an additional demand on the urban w'ater supply. 
Moreover, the addition of trees in the urban setting 
would exert a profound impact on the local 
microclimate and ultimately on the power 
consumption both in summer (as the result of 
reduced air conditioning load) and in winter (due to 
a decreased wind loading on houses). 



These effects are now indirectly quantifiable 
through a combination of remote sensing and 
mathematical modeling of land surface processes. 
We arc currently extending this analysis to examine 
patterns of urbanization in Chester Co, PA, and 
deforestation in Costa Rica. 
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